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We revisit the determination of the pseudo-critical line of QCD with four degenerate quarks at 
non-zero temperature and baryon density by the method of analytic continuation. We determine the 
pseudo-critical couplings at imaginary chemical potentials by high-statistics Monte Carlo simulations 
and reveal deviations from the simple quadratic dependence on the chemical potential visible in 
earlier works on the same subject. Finally, we discuss the implications of our findings for the shape 
of the pseudo-critical line at real chemical potential, comparing different possible extrapolations. 
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I. INTRODUCTION 

The study of QCD at non-zero baryon density by nu- 
merical simulations on a space-time lattice is plagued by 
the well-known sign problem: the fermionic determinant 
is complex and the Monte Carlo sampling becomes un- 
feasible. 

One of the possibilities to circumvent this problem is 
to perform Monte Carlo numerical simulations for imag- 
inary values of the baryon chemical potential, where the 
fermionic determinant is real and the sign problem is ab- 
sent, and to infer the behavior at real chemical potential 
by analytic continuation. 

The idea of formulating a theory at imaginary \i was 
first suggested in Ref. while the effectiveness of the 
method of analytic continuation was pushed forward in 
Ref. @. Since then, the method has been extensively 
applied to QCD % i, i, i, 0, S I, E3 and tested in 
QCD- like theories free of the sign problem (TTI . II2I . lilj 
Qj, GJ, [H, [13 and in spin models pj, GJ]. 

The state-of-the-art is the following: 

• the method is well-founded and works fine within 
the limitations posed by the presence of non- 
analyticities and by the periodicit y o f the theory 
with imaginary chemical potential [20(; 

• the analytic continuation of physical observables is 
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improved if ratios of polynomials (or Pade approx- 
imants [2l[) are used as inter pola ting functions at 
imaginary chemical potential [131 ]; 

• the analytic continuation of the (pseudo-)critical 
line on the temperature - chemical potential plane 
is well- justified, but careful tests in two-color 
QCD 0, [13] and in three-color QCD with finite 
isospin density [I?} have evidenced some difficulties 
in its application; 

• also some partial information about the nature of 
the phase transition as a function of the chemical 
potential can be obtained by a careful study of the 
phase diagram in the T - Im(^) plane @, H2, [Hj]. 

In particular, the numerical analyses in Refs. [3, [ItJ 
have shown that, while there is no doubt that an analytic 
function exists which interpolates numerical data for the 
pseudo-critical couplings for both imaginary and real fi 
across [i = 0, determining this function by an interpola- 
tion of data at imaginary /1 could be misleading. Indeed, 
it was found that non-linear terms in the dependence of 
the pseudocritical coupling /3 C on /i 2 in general cannot be 
neglected and the prediction for the pseudo-critical cou- 
plings at real chemical potentials may be wrong if data at 
imaginary fj, are fitted according to a linear dependence. 
Moreover, the coefficients of the linear and non-linear 
terms in /x 2 in a Taylor expansion of /3 c (/i 2 ) are all nega- 
tive. That often implies subtle cancellations of non-linear 
terms at imaginary chemical potentials (/i 2 < 0) in the 
region available for analytic continuation (first Roberge- 
Weiss sector). The detection of such terms, from simu- 
lations at /i 2 < only, may be difficult and requires an 
extremely high accuracy. 
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In Refs. [TH [ItJ it was realized that, though a poly- 
nomial of order /i 6 seems to be sufficient in all explored 
cases in two-color QCD and in three-color QCD at finite 
isospin density, an increased predictivity can be achieved 
by a fit with the linear term in /i 2 fixed from data at 
small values of fj, 2 only and kept "constrained" when data 
at larger values of /i 2 are included. Moreover, the idea 
was proposed to parameterize the critical line directly in 
physical units in the T, fx plane (instead than in the ft, \x 
plane), and some Ansatze were tested for such parame- 
terization, which provided a very good description of the 
critical line with a reduced number of parameters and an 
increased predictivity. 

The aim of this paper is to apply the experience ac- 
quired through the study of sign-problcm-free theories 
to the determination of the pseudo-critical line in three- 
color QCD with four degenerate quarks. With respect to 
the well-known work of Ref. [IJ , we will take advantage of 
the use of a smaller lattice for collecting more determina- 
tions of pseudo-critical couplings at imaginary chemical 
potential with increased numerical accuracy and will ap- 
ply on the new lattice data the fit strategies worked out 
in our previous studies of Refs. (13. \v\. 



II. NUMERICAL RESULTS 

In our numerical analysis, we consider SU(3) with 
nf = 4 degenerate standard staggered fermions of mass 
am = 0.05 on a 12 3 x 4 lattice. The algorithm adopted 
for Monte Carlo numerical simulations is the exact $ al- 
gorithm described in Ref. {24|, properly modified for the 
inclusion of a finite chemical potential. 
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FIG. 1: Distribution of the real part of the Polyakov loop 
in SU(3) with n/ = 4 on a 12 3 x 4 lattice with am=0.05 at 
a/i = 0.170i and for two ft values around the transition. 

In this theory the critical line in the temperature - 
chemical potential plane is a line of first order transitions, 
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FIG. 2: Monte Carlo history of the real part of the Polyakov 
loop in SU(3) with n/ = 4 on a 12 3 x 4 lattice with am=0.05 
at afi = 0.170i and ,3=5.066. 
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FIG. 3: Susceptibility of the (real part of the) Polyakov loop 
vs ft in SU(3) with nf = 4 on a 12 3 x 4 lattice with am=0.05 
and afi — 0.150i. The solid lines represent the Lorentzian 
interpolation. 



over all the range of fx 2 values in the first Roberge- Weiss 
(RW) sector, -(tt/3) 2 < {^/T) 2 < 0. Tunneling between 
the different phases clearly emerges from the distribu- 
tion on the thermal equilibrium ensemble of the values of 
observables like the (real part of) the Polyakov loop, the 
chiral condensate, the plaquette across the transition (see 
Fig. [T] as an example of the typical two-peak structures). 
As a further evidence of the first order nature of the tran- 
sition, we show in Fig. [2] the Monte Carlo run history of 
the (real part of) the Polyakov loop at afi — 0.1 70i and 
P = 5.066 ~ P C1 which exhibits tunneling events between 
the two phases every few thousands trajectories, on aver- 
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age. Typical statistics have been around 10k trajectories 
of 1 Molecular Dynamics unit for each run, growing up 
to 100k trajectories for 4-5 j3 values around /3 c (/x 2 ), for 
each /i 2 , in order to correctly sample the critical behavior 
at the transition. The critical /3(/j, 2 ) is determined as the 
value for which the susceptibility of the (real part of the) 
Polyakov loop exhibits a peak. To precisely localize the 
peak, a Lorentzian interpolation is used (see Fig. [31 for 
example). In all cases, this kind of determination is com- 
patible with that consisting in estimating the point where 
the two peaks in the distribution of the (real part of the) 
Polyakov loop have equal height, or in locating the peak 
of the susceptibility by the Ferrenberg-Swcndsen method. 
We verified also that the determinations do not change 
if the susceptibility of a different observable, such as the 
baryon number, is used. In Table|T]we summarize our de- 
terminations of the critical couplings: in a few cases we 
have repeated the determination also on a 16 3 x 4 lattice, 
where negligible corrections, within the reported errors, 
have been observed. The plot of the data for /3(/x 2 ) versus 
{afi) 2 - see Fig. |4]- clearly shows that data do not line 
up along a straight line in all the first RW sector, thus 
indicating that the curve /3(/x 2 ) cannot be parametrized 
by a polynomial of order /i 2 . In fact, as we will see soon, 
at least a polynomial of order n 6 is needed to get a fit 
with a reasonable x 2 . 



TABLE I: Summary of the values of /3 c (/x 2 ) for SU(3) with 
n/ = 4 on the 12 3 x 4 lattice with fermionic mass am=0.05. 



a Im(p) 




0. 


5.04259(30) 


0.060 


5.04550(30) 


0.080 


5.04839(30) 


0.100 


5.05121(33) 


0.125 


5.05590(31) 


0.150 


5.06150(30) 


0.170 


5.06647(35) 


0.185 


5.07136(40) 


0.200 


5.07664(30) 


0.210 


5.08031(38) 


0.220 


5.08419(33) 


0.228 


5.08668(30) 


0.235 


5.08961(30) 


0.239 


5.09243(30) 


0.243 


5.09407(30) 


0.2465 


5.09586(30) 


0.250 


5.09754(40) 


0.2535 


5.09970(42) 


0.257 


5.10092(31) 


0.260 


5.10343(30) 


tt/12 


5.1043(5) 



We have tried several kind of interpolations of the crit- 
ical couplings at /i 2 < 0. At first, we have considered 
interpolations with polynomials up to order /x 6 (see Ta- 
ble [II] for a summary of the resulting fit parameters and 
their uncertainties as obtained with the MINUIT min- 
imization code). We can see that data at /x 2 < are 
precise enough to be sensitive to terms beyond the order 
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FIG. 4: Critical couplings obtained in SU(3) with nf = 4 on 
a 12 x 4 lattice with am=0.05. The dashed vertical line in- 
dicates the boundary of the first RW sector, alm(pi) = tt/12. 



/i 2 ; indeed, a good x 2 /d.o.f. is not achieved before in- 
cluding terms up to the order /x 6 . In Fig.[5]Jleft) we show 
how the fit with the 6th-order polynomial compares with 
data of /3 c (/x 2 ); the dotted lines around the fitting curve 
(solid line) delimit the 95% CL band. 

As in Ref. pj}, we performed a "constrained" fit: first, 
the largest interval [(a/x) 2 ^, 0] was identified where data 
could be interpolated by a polynomial in t< 2 , with a 
X 2 /d.o.f ~ 1. It turned out that (a/x) 2 ^ = (0.235x) 2 , 
the quadratic coefficient being equal to —0.8509. Then, 
all available data were fitted by a 6th-order polynomial, 
with the quadratic coefficient fixed at —0.8509 (see Ta- 
ble Hand Fig. fright)). 

Then, we have considered interpolations with ratios of 
polynomials of order up to /x 4 (see Tablc|TT]for a summary 
of the resulting fit parameters). The interpolations with 
the least number of parameters for which we got good fits 
to the data at /x 2 < are the ratio of a 2nd to 4th order 
polynomial and the ratio of a 4th to 2nd order polynomial 
(see Fig. IHJleft) for the latter). 

Finally, we have tried here the fit strategy first sug- 
gested in Ref. [l?} , consisting in writing the interpolating 
function in physical units and to deduce from it the func- 
tional dependence of /3 C on /i 2 , after establishing a suit- 
able correspondence between physical and lattice units. 
The natural, dimensionless variables of our theory are 
T/T c (0), where T c (0) is the critical temperature at zero 
chemical potential, and fi/T. The ratio T/T c (0) is de- 
duced from the relation T — l/(N t a(f3)), where N t is 
the number of lattice sites in the temporal direction and 
a((3) is the lattice spacing at a given /5. Strictly speaking 
the lattice spacing depends also on the bare quark mass, 
which in our runs slightly changes as we change (3 since 
we fix am. However in the following evaluation, which 
is only based on the perturbative 2-loop /3-function, we 
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TABLE II: Parameters of the fits to the critical couplings in SU(3) with n/ = 4 on a 12 3 x 4 lattice with fermionic mass 
am=0.05, according to the fit function /3 C (^ 2 ) = (ao + ai(a/i) 2 + a,2(a^t) 4 + as(a/i) 6 )/(l + a4(ap) 2 + as(a/i) 4 ). Blank columns 
stand for terms not included in the fit. The asterisk denotes a constrained parameter. Fits are performed in the interval 
[(a/i m in) 2 , 0]; the last column gives the value of (a/i m m) 2 . 



Ol 



a 2 
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a 5 x^/d-o-f. (a^min)^ 



6.63 -(tt/12) 2 

0.85 -0.235 2 

2.13 -(tt/12) 2 

1.10 -(vr/12) 2 

1.20 -(tt/12) 2 

1.741(29) 1.13 -(tt/12) 2 

1.09 -(tt/12) 2 



5.04198(22) 
5.04256(24) 
5.04311(36) 
5.04254(50) 
5.04277(27) 
5.04284(28) 
5.04276(27) 



-0.8839(48) 

-0.8509(71) 

-0.761(26) 

-0.892(72) 

-0.8509* 

55.799(14) 

58.196(12) 



1.77(36) 
-3.1(2.4) 
-1.70(55) 

-9.46(13) 



-46.(23.) 
-34.0(8.2) 



11.2266(27) 
11.7044(24) 
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FIG. 5: Fits to the critical couplings: plain 6th order polynomial (left) and 6th-order polynomial with constrained coefficient 
of the quadratic term (right). 
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FIG. 6: Fits to the critical couplings: ratio of a 4th to 2th order polynomial (left) and "physical" fit according to the function @ 
(right). 
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shall neglect such dependence. We use for a(/3) the per- 
turbative 2-loop expression with N c = 3 and rif =4. 
We adopted the 3-parameter function 



Tc(0) 



1 2 



l + C^/T 2 ^) 



l + A^/T^+B^/T^) 



(1) 



leading to the following implicit relation between /3 C and 

a 2 (/3 c (A* 2 ))l2-ioo P = a 2 (/3 c (0))| 2 -ioo P 

l + Af/Tt + B^/T? 

i + <V/t 2 

The values of the fit parameters turned to be 



(2) 



/3 C (0) = 5.04295(25) 
B = 0.12724(75) 



A = 1.00315(95) , 

C = 0.8538(11) , (3) 



with % 2 /d.o.f.=1.26. In Fig. [6](right) we compare the fit 
to data for f3 c (fj, 2 ). 

The important question now is whether the success- 
ful interpolations we found in the fj, 2 < region have a 
consistent extrapolation to fj 2 > 0. In Fig. [7] we have 
plotted the extrapolations to the interval < fi/T < 2 
of the following fits: 

• quadratic fit, performed in the interval — 0.235 2 < 
(an) 2 < (2nd line in Table HJ); 

• sextic constrained polynomial (5th line in Table IPTl); 

• ratio (4,2) of polynomials (last line in Table ITU) ; 

• "physical" fit, Eqs. ©-©. 

The four curves agree as long as fi/T < 0.6, but 
then spread. This means that different interpolations, 
which all reproduce the trend of data in the fit region 
— (7r/12) 2 < (a/i) 2 < and take correctly into account 
the deviation from the quadratic behavior in that region, 
lead to somewhat distinct extrapolations. It is evident 
that, unless an extra-argument is found to make one fit- 
ting function preferable with respect to the others, one 
cannot rely on a unique extrapolation, except in the re- 
gion fi/T < 0.6. The unpleasant aspect is that in the 
same region also deviations from the simple quadratic 
behaviour are negligible: that means that, even if we 
are able to see deviations from the quadratic behaviour, 
we are not able to extrapolate them to real chemical 
potentials in a reliable way, therefore the fact that we 
can see them is in some sense useless. That emerges as 
a shortcoming of analytic continuation, which however 
could be less severe in the more physical case of rif =2 
or rif = 2 + 1. Indeed, in those cases, the curvature of 
the critical line at /i = (i.e. the linear term in fi 2 ) is 
smaller than in the rif = 4 case and larger non-linear 
contributions should be needed to bend the critical line 



towards a critical baryon chemical potential of the or- 
der of 1 GeV at T — 0, so that the sensitivity to such 
non-linear contributions could be hopefully enhanced 1 . 

If one believes that (i) the "physical" fit has a better 
chance to give the correct behavior of the critical line at 
real chemical potential and (ii) it reproduces the critical 
line all the way down to T = 0, then the critical value of 
the chemical potential on the T = axis can be estimated 
as v /Q r B=2.5904(93) T c (0), yielding a critical baryon 
chemical potential at T = slightly above 1 GeV, in 
rough agreement with the expected lightest baryon mass. 
Of course neither of the two assumptions above can be 
supported by valid arguments. 

One may ask whether a different choice of the values of 
the imaginary chemical potential, for which the pseudo- 
critical line has been located, could have improved our 
predictivity. Our choice in the present work has been 
to distribute the values more or less uniformly in fj, 2 : 
we have then increased the density, after a preliminary 
analysis of our data, in the region closer to the RW tran- 
sition, where deviations from the simple linear deviation 
in /i 2 were already visible. Of course the optimal choice 
could be different from that: for instance, once the re- 
gion where the linear approximation works well is known, 
one realizes that fewer points chosen at the border of the 
same region would have provided the same amount of in- 
formation, however this is a hint which is known only a 
posteriori. 

Finally, we present in Fig.[S]an update of Fig. 4 (left) of 
Ref . [26[ , where several determinations of the critical line 
existing in the literature are presented together with the 
results of this work. Looking at Fig. [FJ one could com- 
ment that the extrapolation of the "physical" fit exhibits 
the same trend as data from reweighting, whereas that 
from the sextic constrained fit mimics the strong cou- 
pling behavior 27], the other two extrapolations of ours 
lying in-between. However, previous determinations at 
real fi in the literature seem to be in fair agreement up to 
fi/T ~ 1.2. If one takes this common trend as benchmark 
for our extrapolations, the "physical" and the polynomial 
ratio (4,2) seem to be favoured. 

We have tried to put the previous observation on a 
more solid ground, including in our fit also (subsets of) 
data at real chemical potential available from the litera- 
ture (see Fig. [5J . A serious limitation of this combined 
approach is the inhomogeneity of the data presently avail- 
able, due to different lattices and systematics. Indeed, we 
could not get acceptable values of the x 2 /d.o.f. if not re- 
stricting the region of real chemical potentials included in 
the fit to the interval < a/i ~ 0.6. Unfortunately, this 
is also the region where our extrapolations are consistent 
with each other, so that this combined fit was of little 
help. However, if the inhomogeneity of data at real fi 
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FIG. 7: Extrapolation to real chemical potentials of the quadratic (with (afi) 2 = (0.235i) 2 ), sextic constrained, ratio of 
polynomials and "physical" fits. 
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FIG. 8: Comparison of our extrapolations with other determinations in the literature. For the sake of readability, our extrapo- 
lations have been plotted without error bands and labels, since they can be easily recovered from the previous figure. Legenda: 
D'Elia, Lombardo, Ref. Azcoiti et al, Ref. Fodor, Katz, Ref. [2^|; Kratochvila, de Forcrand, Ref. [2f|. 



will be reduced by new Monte Carlo determinations, the 
combined-fit strategy could bring along an appreciable 
improvement. 



III. DISCUSSION 

In this paper we have revisited the application of the 
method of analytic continuation from imaginary to real 
chemical potential in QCD with rif = 4 degenerate fla- 
vors. The motivations of this study were 

• to determine precisely the pseudo-critical line 
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/3 C (/U 2 ) in the region of negative fi 2 , by sampling 
it through the accurate determination by Monte 
Carlo methods of about 20 data points, almost 
uniformly distributed in the region — (tt/12) 2 < 
(a[i) 2 < and by suitably interpolating these data 
points with an analytic function; 

• to exploit, differently from what was hitherto done 
in the literature, interpolating functions sensitive 
to possible deviations of the critical line from the 
quadratic behavior in /i for larger absolute values of 
fi in the above-mentioned region; these deviations 
were clearly seen in QCD-like theories, such as 2- 
color QCD and finite isospin QCD, where it was 
given compelling evidence that their neglect could 
mislead the analytic continuation to real chemical 
potential; 

• to extrapolate the newly adopted interpolations to 
the region of real /i and to re-determine, therefore, 
the critical line in QCD. 

The outcome has been that deviations from the 
quadratic behavior in /i of the pseudo-critical couplings 
at negative /i 2 are indeed visible in QCD with nf = 4. 
However there are several kinds of functions able to in- 
terpolate them, leading to extrapolations to real /x which 
start diverging from each other for fx/T > 0.6. Unfor- 
tunately in the range of real chemical potentials where 
the different extrapolations agree, deviations from the 
quadratic behavior in fi are negligible, so that our ef- 
forts to determine such deviations in the critical line in a 
consistent way, which were successful for QCD at finite 
isospin density, reveal useless in this case. 

One may ask which significant differences exist be- 
tween QCD at finite isospin density and QCD at finite 
baryon density which could be at the basis of such dif- 
ferent outcomes. Many physical properties distinguish 
the two theories, but the one which is probably most 



relevant to our problem is the different extension of 
the available region for the determination of the criti- 
cal line on the imaginary chemical potential side: for 
QCD with an imaginary baryon chemical potential one 
has < Im(/x)/T < 7r/3, while for QCD with an imagi- 
nary isospin chemical potential one has < Im(/x)/T < 
7r/2 [l7T |. In the second case one has a wider region where 
more information about non-linear terms in /i 2 can be 
collected in order to obtain reliable extrapolations. It 
is indeed interesting to notice that typical deviations of 
/3 c (/z 2 ) — /3 C (0) from the quadratic behaviour, visible at 
the border of the available region on the imaginary chem- 
ical potential side, are of the order of 5% in the present 
work, while they were of the order of 10% for QCD at 
finite isospin chemical potential. 

The shortcomings of the method of analytic continua- 
tion, which emerged in this work on QCD with nf = 4, 
could be less severe in the more physical case of n/ = 2 
or nj = 2 + 1, where the curvature of the critical line 
at /j = (i.e. the linear term in /j, 2 ) is smaller and the 
sensitivity to non-linear terms in /i 2 could be enhanced. 
Substantial improvement could come by either theoret- 
ical developments, able to help discriminating between 
one kind of interpolation and another, or by combined 
numerical strategies, aiming at gathering information on 
the form of the critical line from approaches (such as 
reweighting, Taylor expansion, canonical approach, etc.) 
which have been applied so far only in exclusive way. 
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